#include <math.h>
#include <stdio.h>
#define NRANSI
#include "nrutil.h"

#ifndef TWOPID
#define TWOPID 6.2831853071795865
#endif

float period(double * x, double * y, int n, double ofac, double f_low,
        double f_hi, double * px,double * py, int np, int * nout,
        int * nstart, int * jmax, float * prob)
{
    void avevar(double data[], unsigned long n, float *ave, float *var);

    int i,j;
    float ave,c,cc,cwtau,effm,expy,pnow,pymax,s,ss,sumc,sumcy,sums,sumsh,
    sumsy,swtau,var,wtau,xave,xdif,xmax,xmin,yy;
    double arg,wtemp,*wi,*wpi,*wpr,*wr;
    int percent_done = 0;
    double hifac,f_nyquist,time_interval,lowfac;
    float maxtime,mintime;
    int nout2;

    // some memory assignment ?!?
    wi=dvector(1,n);
    wpi=dvector(1,n);
    wpr=dvector(1,n);
    wr=dvector(1,n);

    // *************************
    maxtime = mintime = x[0];
    for (j = 0; j <= n; j++) {
        if (x[j] < mintime) mintime = x[j];
        if (x[j] > maxtime) maxtime = x[j];
    }

    time_interval = maxtime - mintime;

    f_nyquist = n/(2*time_interval);
    hifac = f_hi/f_nyquist;
    lowfac = f_low/f_nyquist;

    fprintf(stderr,"Time interval: %lf\n",time_interval);
    fprintf(stderr,"Nyquist frequency: %lf\n",f_nyquist);
    // *********************



    // caluculates the number of frequencies in the periodogram
    *nout=0.5*ofac*(hifac)*n;
    nout2=0.5*ofac*(hifac)*n;


    // calculates the required start frequency
    *nstart = 0.5*ofac*lowfac*n;

    fprintf(stderr,"np %d, nout %d, nout2 %d, nstart %d\n\n",
            np,*nout,nout2,*nstart);


    // check to see if the length of output array is sufficient to
    // hold nout frequenies
    if (*nout > np) nrerror("output arrays too short in period");


    // calculates the average and variance
    avevar(y,n,&ave,&var);

    xmax=xmin=x[0];

    fprintf(stderr,"/n** y[0] = %lf\ty[0] = %lf\n",x[0],y[0]);

    for (j=0;j<=n;j++) {
        if (x[j] > xmax) xmax=x[j];
        if (x[j] < xmin) xmin=x[j];
    }

    xdif=xmax-xmin;
    xave=0.5*(xmax+xmin);
    pymax=0.0;
    pnow=1.0/(xdif*ofac);  // frequency


    for (j=0;j<=n;j++) { // looping over time samples
        arg=TWOPID*((x[j]-xave)*pnow); //2pi((t(j)-avetime)*delta_freq)
        wpr[j] = -2.0*SQR(sin(0.5*arg));
        wpi[j]=sin(arg);
        wr[j]=cos(arg);
        wi[j]=wpi[j];
    }

    for (i=0;i<=(*nout);i++) { // looping over frequencies
        px[i]=pnow;   // fills the frequency array
        sumsh=sumc=0.0;
        for (j=0;j<=n;j++) {
            c=wr[j];
            s=wi[j];
            sumsh += s*c;
            sumc += (c-s)*(c+s);
        }
        wtau=0.5*atan2(2.0*sumsh,sumc); // p577 (13.8.5) => w * tau
        swtau=sin(wtau);
        cwtau=cos(wtau);
        sums=sumc=sumsy=sumcy=0.0;

        for (j=0;j<=n;j++)
        {
            s=wi[j];
            c=wr[j];
            ss=s*cwtau-c*swtau;
            cc=c*cwtau+s*swtau;
            sums += ss*ss;
            sumc += cc*cc;
            yy=y[j]-ave;   // (hj-ave_h) => p577
            sumsy += yy*ss;
            sumcy += yy*cc;
            wr[j]=((wtemp=wr[j])*wpr[j]-wi[j]*wpi[j])+wr[j];
            wi[j]=(wi[j]*wpr[j]+wtemp*wpi[j])+wi[j];
        }
        py[i]=0.5*(sumcy*sumcy/sumc+sumsy*sumsy/sums)/var;
        if (py[i] >= pymax) pymax=py[(*jmax=i)];
        pnow += 1.0/(ofac*xdif);
    }
    expy=exp(-pymax);
    effm=(f_hi/f_nyquist)*(*nout)/ofac;
    *prob=effm*expy;
    if (*prob > 0.01) *prob=1.0-pow(1.0-expy,effm);
    free_dvector(wr,1,n);
    free_dvector(wpr,1,n);
    free_dvector(wpi,1,n);
    free_dvector(wi,1,n);
    return(effm);
}

#undef NRANSI
